
;;========================================================================
;; Compare internal distances with each other

mfn = './emaf_paper/save_files/morph_distances.sav'
IF ~ FILE_TEST(mfn) THEN run_morph_matching,/drivemorph,/silent
restore,mfn,/ver

;; Get galactic params
defsysv, '!MW', exists = exists
IF ~ exists THEN galactic_params 
d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART

restore,'./emaf_paper/bgps_emaf.sav',/ver
n = n_elements(emaf)

;;=====================================================================
;; Internal Distance Comparison

myps,'./emaf_paper/plots/internal_dist_comp.eps',xsize=11.5,/cmyk
multiplot_xm,[2,1],xgap=0.030,/xmargin

range = [0, ceil( max(dml[*,0]) > max(dbar[*,0]) )]

plot,[0],[0],/nodata,xr=range,yr=range,/xst,/yst,/iso,$
     xtitle='Maximum Likelihood Distance [kpc]',$
     ytitle='Mean Distance [kpc]'

oploterror,dml[*,0],dbar[*,0],dml[*,1],dbar[*,1],/lobar,$
           errcolor=cgColor('BLK4'),psym=3
oploterror,dml[*,0],dbar[*,0],dml[*,2],dbar[*,1],/hibar,$
           errcolor=cgColor('BLK4'),psym=3

cgplots,range[0],range[0]
cgplots,range[1],range[1],/cont,thick=4,linestyle=2,color='Dodger Blue'

dtan = !MW.R0 * cos(emaf.l * !dtor) / cos(emaf.b * !dtor) / 1.d3
cgloadct,40,/silent,/reverse
color = fix(round((pchoo-0.5)/0.5*255.)<255)>0
psym = 16 - (abs(dml[*,0] - dtan) LE 0.5)*2
symsize = 0.4 + 0.2*(abs(dml[*,0] - dtan) LE 0.5)

for ij=0,n_elements(color)-1 DO $
   cgOplot,dml[ij,0],dbar[ij,0],psym=psym[ij],symsize=symsize[ij],$
           color=color[ij]
           
multiplot,/doyaxis


tan = psym NE 16
near = constr_ml AND ~tan AND dml[*,0] LT dtan
far = constr_ml AND ~tan AND ~near

print,'TAN: ',fix(total(tan))
print,'NEAR:',fix(total(near))
print,'FAR: ',fix(total(far))

tin = where(psym NE 16,ntan)
print,'NTAN: ',ntan
print,'NDOTS: ',n_elements(psym)

cgPlot,abs(dml[*,0]-dbar[*,0]),pchoo,psym=16,symsize=0.4,charsize=1.0,/nodata,$
       xtit='| d!dml!n - d!dbar!n | [kpc]',ytit='P!dchoose!n',yr=[0.45,1.0],/yst

col = replicate('Opposite',n_elements(color))
col[where(psym EQ 14)] = 'Deep Pink'

for ij=0,n_elements(color)-1 DO $
   cgOplot,abs(dml[ij,0]-dbar[ij,0]),pchoo[ij],psym=psym[ij],$
           symsize=symsize[ij],color=col[ij]

vline,/h,0.6827,thick=3,linestyle=3,color='Dodger Blue'

jind = where(pchoo LT 0.6827, njun)
print,'N1sig: ',njun,'  Fraction: ',float(njun)/float(n)
jind = where(pchoo LT 0.9545, njun)
print,'N2sig: ',njun,'  Fraction: ',float(njun)/float(n)
jind = where(pchoo LT 0.8, njun)
print,'N80% : ',njun,'  Fraction: ',float(njun)/float(n)

print,min(pchoo)

cgloadct,40,/silent,/reverse
cgcolorbar,/vert,/right,range=[0.5,1.0],title='P!dchoose!n',charsize=0.8,$
           position=[0.90,!Y.window[0],0.93,!Y.window[1]],divis=0

myps,/done,/mp

;;==============================
;; CDF PLOT

myps,'./emaf_paper/plots/pchoose_cdf.eps'
ph = histogram(pchoo,bin=0.0001,loc=pl)
ph = total(ph,/cum) / total(ph)
cgPlot,pl,ph,xtit='P!dchoose!n',ytit='Cumulative Fraction',thick=3,charsize=1.0
vline,0.6827,thick=3,color='Dodger Blue'
myps,/done


print,m4_stat(psym),m4_stat(symsize)

;;========================================================================
;; Compare distances against the GRS

angsep = 10.

;; Check that the GRS elements of EMAF have been filled in...
IF emaf[0].grs_angsep NE angsep THEN BEGIN
   add_grs_info, ANGSEP=angsep
   undefine,emaf
   restore,'./emaf_paper/bgps_emaf.sav',/ver
   n = n_elements(emaf)
ENDIF  

;; Restrict consideration to "well-constrained" sources...
cind  = where(constr_ml)
duse   = duse[cind,*]
dbar  = dbar[cind,*]
emaf  = emaf[cind]
pchoo = pchoo[cind]
color = color[cind]


;; SELECT ONLY GRS-ASSOCIATED SOURCES
ind = where(emaf.grs, n)
emaf = emaf[ind]
duse = duse[ind,*]
dbar = dbar[ind,*]
pchoo = pchoo[ind]
color = color[ind]


;;==================================================
;; MAKE PLOTS!!1!
myps,'./emaf_paper/plots/grs_dist_comp.eps',xsize=10,/cmyk
multiplot_xm,[2,1],/xmargin

plot,[0],[0],/nodata,xr=range,yr=range,/xst,/yst,/iso,$
     xtitle='Maximum Likelihood Distance [kpc]',$
     ytitle='GRS Distance [kpc]'

oploterror,duse[*,0],emaf.grs_dist,duse[*,1],fltarr(n),/lobar,$
           errcolor=cgColor('BLK4'),psym=3
oploterror,duse[*,0],emaf.grs_dist,duse[*,2],fltarr(n),/hibar,$
           errcolor=cgColor('BLK4'),psym=3

cgplots,range[0],range[0]
cgplots,range[1],range[1],/cont,thick=4,linestyle=2,color='TG4'

cgloadct,40,/silent,/reverse
;color = fix(round(emaf.c/0.3*255.))<255
FOR ii=0,n-1 DO $
   cgPlots,duse[ii,0],emaf[ii].grs_dist,psym=16,symsize=0.6,$
         color=color[ii]

cgColorbar,range=[0.5,1.0],/yst,/vert,/right,tit='P!dchoose!n',charsize=0.9,$
           pos=[0.89,0.17,0.92,0.90],divis=0

bind = where( abs(duse[*,0]-emaf.grs_dist) GT 1., nbd)
print,'Nbad:',nbd,' Fraction Matching: ',1.d - float(nbd)/float(n)

multiplot

plot,[0],[0],/nodata,xr=range,yr=range,/xst,/yst,/iso,$
     xtitle='Mean Distance [kpc]'

oploterror,dbar[*,0],emaf.grs_dist,dbar[*,1],fltarr(n),/lobar,$
           errcolor=cgColor('BLK4'),psym=3
oploterror,dbar[*,0],emaf.grs_dist,dbar[*,1],fltarr(n),/hibar,$
           errcolor=cgColor('BLK4'),psym=3

cgplots,range[0],range[0]
cgplots,range[1],range[1],/cont,thick=4,linestyle=2,color='TG4'

cgloadct,40,/silent,/reverse
FOR ii=0,n-1 DO $
   cgPlots,dbar[ii,0],emaf[ii].grs_dist,psym=16,symsize=0.6,$
         color=color[ii]


myps,/done,/mp

;;==================================================================
;; Make another one, but this time for "/ dtan"

tandist = !MW.R0*cos(emaf.l * !dtor)/cos(emaf.b * !dtor)/1.d3

FOR j=0,2 DO duse[*,j] /= tandist
FOR j=0,1 DO dbar[*,j] /= tandist
emaf.grs_dist /= tandist
range=[0,2]

;; myps,'./emaf_paper/plots/grs_dist_comp2.eps',xsize=10,/cmyk
myps,'./emaf_paper/plots/grs_dist_comp2_pchoo.eps',xsize=10,/cmyk
multiplot_xm,[2,1],/xmargin

plot,[0],[0],/nodata,xr=[range[0],range[1]-0.00001],yr=range,/xst,/yst,/iso,$
     xtitle='ML Distance / d!dtan!n',$
     ytitle='GRS Distance / d!dtan!n'

oploterror,duse[*,0],emaf.grs_dist,duse[*,1],fltarr(n),/lobar,$
           errcolor=cgColor('BLK4'),psym=3
oploterror,duse[*,0],emaf.grs_dist,duse[*,2],fltarr(n),/hibar,$
           errcolor=cgColor('BLK4'),psym=3

cgplots,range[0],range[0]
cgplots,range[1],range[1],/cont,thick=4,linestyle=2,color='TG4'

cgloadct,40,/silent,/reverse
;; color = fix(round( (emaf.l - 15)/50.25*254.)< 254)>0
color = fix(round( (pchoo-0.5)/0.5*254.)< 254)>0
;; color = fix(round( abs(emaf.b)/0.5*254.)< 254)>0

FOR ii=0,n-1 DO $
   cgPlots,duse[ii,0],emaf[ii].grs_dist,psym=16,symsize=0.6,$
         color=color[ii]

;; cgColorbar,range=[0,0.5],/vert,/right,tit='abs(GLAT) [deg]',charsize=0.9,$
cgColorbar,range=[0.5,1.0],/vert,/right,tit='P!dchoose!n',charsize=0.9,$
           pos=[0.89,!Y.window[0],0.92,!Y.window[1]],divis=0


bind = where( abs(duse[*,0]-emaf[*].grs_dist) GT 0.2 AND pchoo GE 0.90, nbd)
print,nbd
for jji=0,nbd-1 do $
   print,emaf[bind[jji]].cnum,color[bind[jji]]
psym = psym * 0. + 16
psym[bind] = 14


multiplot

plot,[0],[0],/nodata,xr=range,yr=range,/xst,/yst,/iso,$
     xtitle='Mean Distance / d!dtan!n'

oploterror,dbar[*,0],emaf.grs_dist,dbar[*,1],fltarr(n),/lobar,$
           errcolor=cgColor('BLK4'),psym=3
oploterror,dbar[*,0],emaf.grs_dist,dbar[*,1],fltarr(n),/hibar,$
           errcolor=cgColor('BLK4'),psym=3

cgplots,range[0],range[0]
cgplots,range[1],range[1],/cont,thick=4,linestyle=2,color='TG4'

cgloadct,40,/silent,/reverse
FOR ii=0,n-1 DO $
   cgPlots,dbar[ii,0],emaf[ii].grs_dist,psym=psym[ii],symsize=0.6,$
           color=color[ii]

myps,/done,/mp



END
